!-------------------------------------------------------------------------------------
subroutine evolve_F_collisions(dfu, dt)
!-------------------------------------------------------------------------------------

  use dimensions

  implicit none

  real,    intent(inout) :: dfu(jv, jz, jy, jx, jdir)
  real,    intent(in)    :: dt


  real                   ::  s(jv, jz, jy, jx, jdir)
  real                   :: s0(jv, jz, jy, jx, jdir)
  real                   :: df(jv, jz, jy, jx, jdir)
 

  df = dfu;

  call rk_rhs_collisions(df, s)
  df = dfu + s*dt/2
  s0 = s

  call rk_rhs_collisions(df, s)
  df = dfu + s*dt/2
  s0 = s0 + 2*s

  call rk_rhs_collisions(df, s)
  df = dfu + s*dt
  s0 = s0 + 2*s

  call rk_rhs_collisions(df, s)
  s0 = s0 + s

  dfu = dfu + s0*dt/6

end subroutine evolve_F_collisions


!-------------------------------------------------------------------------------------
subroutine rk_rhs_collisions(f, rhs)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, vPhase, CollisionRate
  use fft_wrap

  implicit none

  real,    intent(in)    :: f(jv, jz, jy, jx, jdir)
  real,    intent(out)   :: rhs(jv, jz, jy, jx, jdir)
 
  integer                :: ix, iy, iz, iDir
  complex                :: fv(jv)

  if (CollisionRate .eq. 0) return;

  do iDir=1,jDir
     do ix=1, jx
        do iy=1, jy
           do iz=1, jz
           
              fv(:) = f(:, iz, iy, ix, iDir) * (vvec + vPhase)

              call fft_forward(jv,fv,fv)
              fv = fv/float(jv)
              fv(:) = fv(:) * iunit * kVvec(:) 
              call fft_backward(jv,fv,fv)

              rhs(:, iz, iy, ix, iDir) = fv(:) * CollisionRate
           end do
        end do
     enddo
  enddo

end subroutine rk_rhs_collisions

!-------------------------------------------------------------------------------------
